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1 ABSTRACT 

This paper advances the state-of-the-art in 
spray computations with some of our recent contri- 
butions involving scalar Monte Carlo PDF (Proba- 
bility Density Function), unstructured grids and par- 
allel computing. It provides a complete overview of 
the scalar Monte Carlo PDF and Lagrangian spray 
computer codes developed for application with un- 
structured grids and parallel computing. Detailed 
comparisons for the case of a reacting non-swirling 
spray clearly highlight the important role that chem- 
istry/turbulence interactions play in the modeling of 
reacting sprays. The results from the PDF and non- 
PDF methods were found to be markedly different 
and the PDF solution is closer to the reported exper- 
imental data. The PDF computations predict that 
some of the combustion occurs in a predominantly 
premixed-flame environment and the rest in a pre- 
dominantly diffusion-flame environment. However, 
the non-PDF solution predicts wrongly for the com- 
bustion to occur in a vaporization-controlled regime. 
Near the premixed flame, the Monte Carlo particle 
temperature distribution shows two distinct peaks: 
one centered around the flame temperature and the 
other around the surrounding-gas temperature. Near 
the diffusion flame, the Monte Carlo particle tem- 
perature distribution shows a single peak. In both 
cases, the computed PDF’s shape and strength are 
found to vary substantially depending upon the prox- 
imity to the flame surface. The results bring to the 
fore some of the deficiencies associated with the use 
of assumed-shape PDF methods in spray computa- 
tions. Finally, we end the paper by demonstrating 
the computational viability of the present solution 
procedure for its use in 3D combustor calculations 
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by summarizing the results of a 3D test case with 
periodic boundary conditions. For the 3D case, the 
parallel performance of all the three solvers (CFD, 
PDF, and spray) has been found to be good when 
the computations were performed on a 24-processor 
SGI Origin work-station. 

2 NOMENCLATURE 

A pre-exponent of an Arrhenius 

reaction rate term 

a non-unity exponent of an Arrh- 

enius reaction rate term 
a n outward area normal vector 

of the nth surface, m 2 
Bk Spalding transfer number 

b non-unity exponent of an 

Arrhenius reaction rate term 
Cd drag coefficient 

C p specific heat, J/(kg K) 

C<p a constant in Eq. (39) 

c n convection/diffusion coefficient 

of the nth face, kg/s 

D turbulent diffusion coefficient, m 2 /s 

d droplet diameter, m 

E a activation energy of an 

Arrhenius reaction rate term 
h specific enthalpy, J/kg 

Jf diffusive mass flux vector, kg/ms 

k turbulence kinetic energy, m 2 /s 2 

Ik latent heat of evaporation, J/kg 

h,eff effective latent heat of evaporation, 

J/kg (defined in Eq. (11)) 

Mi molecular weight of ith 

species, kg/kg-mole 
mk droplet vaporization rate, kg/s 

rriko initial mass flow rate associated 

with kth droplet group, kg/s 
N av number of time steps employed in 

the PDF time-averaging scheme 
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Nf number of surfaces contained in 
a given computational cell 
N m total number of Monte Carlo 
particles per grid cell 

N p total number of computational cells 
Tik number of droplets in kth group 
P pressure, N/m 2 

P T Prandtl number 

p joint scalar PDF 

Ru gas constant, J / (kg K) 

Re Reynolds number 

r* droplet radius, m 

rko initial droplet radial location, m 
Sk droplet radius squared, r 2 , m 2 
$mic liquid source contribution of the 
gas-phase continuity equation 
Smie liquid source contribution of the 
gas-phase energy equation 
5 m / m liquid source contribution of the 
gas-phase momentum equations 
s m \ s liquid source contribution of the 
gas-phase species equations 
s a liquid source contribution of the a variable 
T temperature, K 

t time, s 

Ui ith velocity component, m/s 

Uik ith velocity component of 

kth droplet group, m/s 
V c volume of the computational cell, m 3 
w a chemical reaction rate, 1/s 

Wj gas-phase chemical reaction rate, 1/s 

Xi Cartesian coordinate in the ith direction, m 
yj mass fraction of jth species 

x spatial vector 

X , mole fraction 

At local time step used in the 

PDF computations, s 

A tj local time step in the flow solver, s 

A tg\ global time step in the spray solver, s 

A tu fuel injection time step, s 

At m j allowable time step in the spray solver, s 
AV computational cell volume, m 3 

6 Dirac-delta function 

e rate of turbulence dissipation, m 2 /s 3 

€j species mass fraction at the droplet surface 

ects species mass fraction at the droplet surface 

r* turbulent diffusion coefficient, kg/ms 

A thermal conductivity, J/(ms K) 

p dynamic viscosity, kg/ms 

u> turbulence frequency, 1/s 


<j> represents a set of scalars of the joint PDF 
independent composition space 
p density, kg/m 3 
a dimensionality of ^-space 
r stress tensor term, kg/ms 2 
9 void fraction 

Subscripts 

/ represents conditions associated with fuel 
g global or gas-phase 

i index for the coordinate or species components 
j index for the species component 
k droplet group or liquid-phase 
/ liquid-phase 

m conditions associated with N m 
n nth-face of the computational cell 
o initial conditions or oxidizer 
p conditions associated with the 
properties of a grid cell 
s represents conditions at the droplet 
surface or adjacent computational cell 
t conditions associated with time 
a index for the scalar component of 
the joint PDF equation 
, partial differentiation with respect 
to the variable followed by it 

Superscripts 

Favre averaging 

time averaging or average based on the Monte 
Carlo particles present in a given cell 
// fluctuations 

3 INTRODUCTION 

Spray combustion is of interest in a wide- 
variety of applications: gas-turbine combustors, in- 
ternal combustion engines (diesel and spark-ignition), 
liquid-rocket motors, and industrial furnaces. A large 
fraction of the world’s energy needs are met by the 
combustion of the liquid fuels. Because of its abun- 
dant use, the need for developing better predictive 
tools for aiding in the design of ever more efficient, 
pollution-free, and stable combustion devices has re- 
ceived considerable attention. 1 ” 10 

The physical modeling of turbulent spray flames 
requires consideration of various complex and rate- 
controlling processes associated with turbulent trans- 
port, mixing, and chemical kinetics, fluid-dynamic 
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characteristics of fuel injection and spray formation, 
transport and vaporization characteristics of individ- 
ual droplets, and the interaction of turbulence with 
chemical kinetics, among others. The interaction of 
turbulence with chemical kinetics may occur over 
a wide-range of disparate time and length scales. 
Turbulence plays an important role in determining 
the rates of mass and heat transfer, chemical reac- 
tions, and liquid-phase evaporation in many practi- 
cal combustion devices. The influence of turbulence 
in a gaseous diffusion flame manifests itself in sev- 
eral forms, ranging from the so-called wrinkled or 
stretched flamelets regime to the distributed combus- 
tion regime, depending upon how turbulence inter- 
acts with various flame scales. 11-12 Although most 
of the turbulent spray combustion models are based 
on either diffusion or premixed gaseous flame theo- 
ries, combustion in a spray flame is more complex 
and seldom occurs in a single mode. 

There are several approaches used in the study 
of turbulent spray flames. 1-10 El Banhawy and 
Whitelaw 2 in their study on the prediction of the 
local flow properties in a spray flame employed an 
assumed shape PDF model with values for its mean 
and variance obtained from the solution of additional 
transport equations. Raju and Sirignano 4 in their 
study on multicomponent spray computations in a 
swirl-stabilized center-body combustor employed the 
eddy break-up model of Spalding 13 to take some ef- 
fect of turbulence on combustion. The eddy break-up 
model might be applicable in cases where chemical 
reaction rates are either very fast or slow compared 
with the turbulent time scales. 

Most of the turbulence closure models for reac- 
tive flows have difficulty in treating nonlinear reaction 
rates. 11-12 The use of assumed shape PDF methods 
was found to provide reasonable predictions for pat- 
tern factors and NO* emissions at the combustor 
exit. 14 However, their extension to multi-scalar chem- 
istry becomes quite intractable. The solution proce- 
dure based on the modeled joint composition PDF 
transport equation has an advantage in that it treats 
the nonlinear reaction rates without any approxima- 
tion. This approach holds the promise of modeling 
various important combustion phenomena relevant 
to practical combustion devices such as flame extinc- 
tion and blow-off limits, and unburnt hydrocarbons 
(UHC), CO, and NOx predictions. 11-12,14 

The success of any numerical tools used in the 
multidimensional combustor modeling depends not 
only on the modeling and numerical accuracy con- 


siderations but also depends on the computational 
efficiency considerations as determined by the com- 
puter memory and turnaround times afforded by the 
present-day computers. With the aim of developing 
an efficient solution procedure for use in multidimen- 
sional combustor modeling, we extended the scalar 
Monte Carlo PDF approach to the modeling of sprays 
with parallel computing in order to facilitate large- 
scale combustor applications. 7 In this approach, the 
mean gas-phase velocity and turbulence fields are de- 
termined from the solution of a conventional CFD 
method, the scalar fields of species and enthalpy from 
a modeled PDF transport equation using a Monte 
Carlo method, and a Lagrangian-based dilute spray 
model is used for the liquid-phase representation. 
The application of this method showed reasonable 
agreement when detailed comparisons were made for 
two different cases involving an open and a confined 
swirl-stabilized spray flames. 7 

It is well known that considerable effort usually 
goes into generating structured-grid meshes for grid- 
ding up practical combustor geometries which tend 
to be very complex in shape and configuration. The 
grid generation time could be reduced considerably 
by making use of existing, automated unstructured 
grid generators. 15-18 With the aim of advancing the 
current multi-dimensional computational tools used 
in the design of advanced technology combustors, two 
new computer codes - LSPRAY 8 and EUPDF 19 - 
were developed here, thereby extending our previ- 
ous work 7 on the Monte Carlo PDF and sprays to 
unstructured grids as a part of the National Com- 
bustion Code (NCC) activity. NCC is being devel- 
oped in the form of a collaborative effort between 
NASA GRC, aircraft engine manufacturers, and sev- 
eral other government agencies and contractors. 17 
The unstructured 3D solver is designed to be mas- 
sively parallel and accommodates the use of an un- 
structured mesh with mixed elements comprising of 
either triangular, quadrilateral, and/or tetrahedral 
type. The ability to perform the computations on 
unstructured meshes allows representation of com- 
plex geometries with relative ease. The application of 
the unstructured grid extension to a confined swirl- 
stabilized spray flame provided reasonable agreement 
with the available droplet velocity measurements. 9 

A current status of the the use of the paral- 
lel computing in turbulent reacting flows involving 
sprays, scalar Monte Carlo PDF and unstructured 
grids was described in Ref. 10. It also outlines sev- 
eral numerical techniques developed for overcoming 
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some of the high computer time-and-storage limita- 
tions placed by the use of Monte Carlo solution meth- 
ods. The parallel performance of both the PDF and 
CFD computations was found to be excellent but the 
results were mixed for the spray module showing rea- 
sonable performance on massively parallel computers 
like Cray T3D; but its performance was poor on the 
workstation clusters. In order to improve the par- 
allel performance of the spray module, two different 
domain decomposition strategies were developed and 
the results from both strategies were summarized. 10 

The main objective of our present work is to 
to investigate the importance of considering chem- 
istry/turbulence interactions in the calculation of 
a reacting spray. This was done by making de- 
tailed comparisons for the case of a reacting non- 
swirling spray for which experimental data was re- 
ported by McDonell and Samuelson. 20 The compar- 
isons involved predictions from two different sets of 
computations, one in which the solution for the tem- 
perature and species fields is obtained from the use of 
the scalar Monte Carlo PDF method and in the other 
they are obtained from the solution of a conventional 
CFD solution. The second objective is to demon- 
strate the computational viability of the present so- 
lution procedure for its use in 3D computations. This 
was done by summarizing the results of a 3D test case 
which was designed primarily to examine the appli- 
cability of the spray particle search algorithm in 3D 
computations and of the newly-implemented periodic 
boundary conditions. For clarity, the overall solution 
procedure from Ref. 10 is repeated here. However, 
for a detailed account of the parallel performance to- 
gether with the development and implementation of 
the parallel method, the interested reader is referred 
to Ref. 10. 

First, complete details of the overall solution 
procedure with a particular emphasis on the PDF 
and spray algorithms are presented along with sev- 
eral other numerical issues related to the coupling 
between the CFD, spray, and PDF solvers. It is fol- 
lowed by the results-and-discussion section where the 
application of the method to predict the local prop- 
erties of two different cases are presented. Finally, we 
conclude the paper by presenting a brief summary of 
important results along with the parallel performance 
of the three solvers (CFD, PDF, and spray). 

4 GOVERNING EQUATIONS FOR THE 
GAS PHASE 

Here, we summarize the conservation equations for 


the gas-phase in Eulerian coordinates derived for the 
multicontinua approach. 21 This is done for the pur- 
pose of identifying the interphase source terms arising 
from the exchanges of mass, momentum, and energy 
with the liquid-phase. 

The conservation of the mass leads to: 

[pVch + \pV c Ui\ tXi = s m ic n k m k ( 1 ) 

k 

For the conservation of the jth species, we have: 

fflcVjh + h oV c Uiyj],x, - [ pV c Dy jiXi ] iXi - pV c wj = 

Smis = €j n k m k (2) 

k 

where 

=0anc2 = 1 

j i 

For the momentum conservation, we have: 

\pVcUi],t + \pVcUiUj] }Xi + \pV e \ g- 

j\,Xj ~ [(1 ^)^c T lii\,Xj = Smlm — 

n k m k u ki -V ^ Pkr\n k u ki[t (3) 

* k * 

where 0 = the void fraction of the gas which is defined 
as the ratio of the equivalent volume of gas to a given 
volume of a gas and liquid mixture. For dilute sprays, 
the void fraction is assumed to be equal to one. The 
shear stress r, ij in Eq. (3) is given by: 

T ij = + Uj^Xi] “ 0 

For the energy conservation, we have: 

[pV c hl t +[pV c Uihl x -[(l-0)V c A f T r J,^ 

— [$VcP],t = Sm/e = 'y ^ 77ljb — lk f ejf' S j (4) 
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5 GAS-PHASE SCALAR JOINT PDF 
EQUATION 

The transport equation for the density-weighted joint 
PDF of the compositions, p, is: 

Wi.t + [?«*£],*, 

{ Transient } {Mean convection } 

+ [pw a (f)p\ tiPa = -[p < u” I rp > p\ )Xl 
{Chemical reactions} {Turbulent convection } 

~[p< ~ p Ji, t, I ± > ft ,* . -\p<y a \±> p],^ a (5) 
{Molecular mixing } {Liquid— phase contribution} 
where 

w a — chemical source term 

for the o-th 
composition variable, 

< u i I = conditional average of 

Favre velocity 
fluctuations, 

< I $ > = conditional average 

of scalar dissipation, 
and 

< j > = conditional average of 

spray source terms. 

The terms on the left-hand side of the above 
equation could be evaluated without any approxima- 
tion, but the terms on the right-hand side of the equa- 
tion require modeling. The first term on the right 
represents transport in physical space due to turbu- 
lent convection. 12 Since the joint PDF, p, contains 
no information on velocity, the conditional expecta- 
tion of < u” | > needs to be modeled. It is mod- 

eled based on a gradient-diffusion model with infor- 
mation supplied on the turbulent flow field from the 
flow solver. 12 

- < u" \ ±> p = T tP't, (6) 

The fact that the turbulent convection is modeled 
as a gradient-diffusion makes the turbulent model no 
better than the k — e model. The uncertainties asso- 
ciated the use of a standard k — c turbulence model to 
swirling flows are well known. 22 Some of the modeling 
uncertainties associated with the use of the standard 
k — e model would be addressed in our future studies 
with the implementation of a non-linear k - e devel- 
oped for the modeling of swirling flows. 22 


The second term on the right-hand side repre- 
sents transport in the scalar space due to molecu- 
lar mixing. A mathematical description of the mix- 
ing process is rather complicated, and the interested 
reader is referred to Ref. 12. Molecular mixing is 
accounted for by making use of the relaxation to the 
ensemble mean submodel. 11 


< -J? x . I ^ >= -C>(^ - $*) (7) 

P “ 

where u> = e/&, and C$ is a constant. For a 

conserved scalar in a homogeneous turbulence, this 
model preserves the PDF shape during its decay, but 
there is no relaxation to a Gaussian distribution. 12 
However, the results of Ref. 14 indicate that the 
choice between the different widely-used mixing mod- 
els is not critical in the distributed reaction regime of 
premixed combustion as long as the turbulent mixing 
frequencies are above 1000 Hz. Most of the practical 
combustors seem to operate at in-flame mixing fre- 
quencies of 1000 Hz and above. The application of 
this mixing model seemed to provide some satisfac- 
tory results when applied to flows representative of 
those encountered in the gas-turbine combustion. 14 

The third term on the right-hand side represents 
the contribution from the spray source terms: 

< ~ p S<* I ± >= ^ nk mk ( 8 ) 

where <f> a = y a , a = 1, 2 , s — a — l 



C pi = Zz(Au + A 2i T + A Zi T 2 + A 4 ;T 3 + A 5i T A ), 

hjt is the heat of formation of ith species, Ru is the 
universal gas constant, e as is a mass fraction of the 
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evaporating species at the droplet surface, and h } eff 
is the effective latent heat of vaporization as modified 
by the heat loss to the droplet interior: 


(15) 


C D = 


24 

Rek 
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l k J}= I k+ 4ir^ 

m k 



( 11 ) 


Here we assumed that the spray source terms 
could be evaluated independent of the fluctuations in 
the gas-phase compositions of species and enthalpy. 
Eqs. (8)-(10) represent the modeled representation 
for the conditional averages of the spray contribution 
to the PDF transport equation. 


6 LIQUID-PHASE EQUATIONS 


For droplet size, the droplet regression rate is 
determined from one of three different correlations 
depending upon the droplet-Reynolds- number range. 
When Rek > 20, the regression rate is determined 
based on a gas-phase boundary-layer analysis valid 
for Reynolds numbers in the intermediate range. 23 
The other two correlations, valid when Re k < 20, are 
taken from Clift et al. 24 


dsh_ 

dt 


- 2 — 
Pk L 


2 „ 

-Re k 

7T 


-.1/2 


f(Bk) 


if Re k > 20 


The spray model is based on the multicontinua ap- 
proach which allows for resolution on a scale greater 
than the average spacing between two neighboring 
droplets. A Lagrangian scheme is used for the liquid- 
phase equations as it eliminates errors associated 
with numerical diffusion. The vaporization model of a 
polydisperse spray takes into account the transient ef- 
fects associated with the droplet internal heating, the 
forced convection effects associated with droplet in- 
ternal circulation and the phenomena associated with 
boundary layers and wakes formed in the intermedi- 
ate droplet Reynolds number range. 4 The present for- 
mulation is based on a deterministic particle-tracking 
method and on a dilute spray approximation which 
is applicable for flows where the droplet loading is 
low. Not considered in the present formulation are 
the effects associated with the droplet breakup, the 
droplet/shock interaction, the multi-component na- 
ture of liquid spray and the phenomena associated 
with dense spray effects and super-critical conditions. 
The spray method provided some favorable results 
when applied to several unsteady and steady-state 
calculations. 4-10 

For the particle position of the kth droplet 
group, we have: 


dxik 

~df 


= u ik 


( 12 ) 


^ [l + (1 + fle*) 1/3 ] Re° k 077 ln( 1 + B k ) 

if 1 <Re k < 20 (16) 

^ ~ [l + (1 + Re k ) 1/3 ] ln(l+B k ) if Re k < 1 

where B k is the Spalding transfer number defined in 
Eq. (22). The function f(B k ) is obtained from the 
solution of Emmon’s problem. 25 The range of validity 
of this function was extended in Raju and Sirignano 4 
to consider the effects of droplet condensation. 

The internal droplet temperature is determined 
based on a vortex model. 23 The governing equation 
for the internal droplet temperature is given by: 


dt 


= 17 


A* 


CpiPirl L da 2 


d 2 T k 


a- 


+ ( 1 + C(<)»)^- 


(17) 


where 



where a represents the coordinate normal to the 
streamsurface of a Hill’s Vortex in the circulating 
fluid, and C(t) represents a nondimensional form of 
the droplet regression rate. The initial and boundary 
conditions for Eq. (17) are given by: 


For the droplet velocity: 

d _ 3 Cppgsjlek 


dt 


16 


Pkr\ 


[Uig ~ U.jfe] 


where 


Re k =2 r -^[(u g -u k ).(u g -u k )f 7 

Pgs 


t — tinjectioni T k — X&, 


(13) 

P 

II 

o 

(14) 

f^l <3 

ii 

<3 


m 

da 


17 


Cptpi 


A, 


> m 

: dt 


3 p k 
'32 A, 


<MEuzIbil h 


B k 


(19) 

( 20 ) 

dsk 

dt 

( 21 ) 
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where a = 0 refers to the vortex center, and a = 1 
refers to the droplet surface. 

The Spalding transfer number is given by: 

R _ Cp( T 9 - _ (yjj, - yj) 

k ~ h.efj (i -yj>) 1 j 

yi ' = 1 + w, (*/'• - 0 < 23) 

where Af 0 is the molecular weight of the gas excluding 
fuel vapor. 

Based on the assumption that phase equilibrium 
exists at the droplet surface, the Clausius- Clap eyron 
relationship yields 



Ik 

(1 

JA1 

,R U 

u 

T k ,)_ 


(24) 


In Eq. (14), the molecular viscosity is evaluated at a 
reference temperature using Sutherland’s equation 


p„(T r ' J ) = 1.4637 10 - 6 


rf 3/2 
_ ref 


T re f + 120 


(25) 


where 

Tref = \T g + | T k) (26) 

The droplets may evaporate, move along the 
wall surfaces, and/or reflect with reduced momen- 
tum upon droplet impingement with the combustor 
walls. In our present computations, subsequent to 
the droplet impingement with the walls, the droplets 
are assumed to flow along the wall surfaces with a 
velocity equal to that of the surrounding gas. 

7 DETAILS OF DROPLET FUEL 
INJECTION 

The success of any spray model depends a great 
deal on the specification of the appropriate injector 
exit conditions. However, a discussion involving the 
physics of liquid atomization is beyond the scope of 
this subject matter. In our present computations, the 
liquid fuel injection is simulated by introducing a dis- 
cretized parcel of liquid mass in the form of spherical 
droplets at the beginning of every fuel-injection time 
step. 

For certain cases, the fuel-injection time step, 
A tip needs to be determined based on the resolution 
permitted by the length and time scales associated 
with several governing parameters such as average 


grid spacing and average droplet spacing and velocity. 
However, our experience showed that for the case of a 
steady-state solution, a time step based on the aver- 
age droplet lifetime yields better convergence . 4 “ 6 Its 
value typically ranges between 1 and 2 milli-seconds 
for the case of reacting flows. 

The spray computations facilitate fuel injection 
through the use of a single fuel injector comprised of 
different holes . 5-6 However, multiple fuel injection in 
a steady-state calculation could be simulated by sim- 
ply assigning different initial conditions for the spatial 
locations of the droplet groups associated with each 
one of the different holes. For a polydisperse spray, 
the spray computations require inputs for the num- 
ber of droplet groups in a given stream and for the 
initial droplet locations and velocities. However, the 
number of droplets in a given group and their sizes 
could be either input directly or computed from a 
properly chosen function for the droplet size distribu- 
tion. The specified initial inputs should be represen- 
tative of the integrated averages of the experimental 
conditions . 5-10 

One correlation typical of those used for the 
droplet size distribution is taken from Ref. 2: 


— = 4.21 10 6 
n 



16.98(4) 


dd_ 

^32 


(27) 


where n is the total number of droplets and dn is 
the number of droplets in the size range between d 
and d + dd. The Sauter mean diameter, ^ 32 , could 
be either specified or estimated from the following 
correlation : 26 


C?32 = Bd 


2ncri 

774 


AI 


(28) 


where Bd is a constant, Vt is the average relative 
velocity between the liquid interface and the ambi- 
ent gas, and A£, is a function of the Taylor number, 

{pl^f)/(P S pf4)- 

A typical droplet size distribution obtained 
from the above correlation in terms of the cumulative 
percentage of droplet number and mass as a function 
of the droplet diameter is shown in Fig. I . 7 
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Figure 1 Droplet-size distribution. 

8 CFD SOLUTION ALGORITHM 

The gas-phase mass and momentum conservation 
equations together with the standard k — c turbulence 
equations with wall functions are solved by making 
of a finite- volume solver with an explicit fourth-stage 
Runge-Kutta scheme. Further details of the code can 
be found in Refs. 16 & 18. 

9 PDF SOLUTION ALGORITHM 

In order to facilitate the integration of the Monte 
Carlo PDF method in a finite-volume context, the 
volume integrals of convection and diffusion in Eq. 
(5) were first recast into surface integrals by means of 
a Gauss’s theorem. 27 Partial integration of the PDF 
transport equation would yield: 

c A i 

P p (£,t + At) = (1 - J£y)Pp{±,t) 

- At i w c(±)p\,i, a 

- A <[< ~ p Ji,x, I ± > ~ A *[< | rp_ > p\^ a 

(29) 


where subscript n refers to the nth-face of the com- 
putational cell. The coefficient c n represents the 
transport by convection and diffusion through the 
nth-face of the computational cell, p . The convec- 
tion/diffusion coefficients in the above equation are 
determined by one of the following two expressions: 

=r^( A ^^^ ) + maz[0,-pa n .M n ] 

2 a CL 

c n = max[\Q.hpa n .uj , I»( ) ] - 0.5 p^.j^ 
and 

c p ~ 5 -/ Cn 

n 

In both the above expressions for c n , a cell-centered 
finite-volume derivative is used to describe the vis- 
cous fluxes; but an upwind differencing scheme is used 
for the convective fluxes in the first expression and a 
hybrid differencing scheme in the second. 

9.1 Numerical Method Based on 
Approximate Factorization 

The transport equation is solved by making use of an 
approximate factorization scheme. 12 Eq. (29) can be 
recast as: 

P P {±, t + At) = 

(I+AtR)((I+AtS)(I+AtM)(I+AtT)p p t)+0(At 2 ) 

(30) 

where I represents the unity operator and T, M, S, 
and R denote the operators associated with spatial 
transport, molecular mixing, spray, and chemical re- 
actions, respectively. The operator is further split 
into a sequence of intermediate steps: 

%(±,t) = {I + AtT)p F {±,t) (31) 

p“{±,t) = {I+AtM)f p {±,t) (32) 

P?*(±, <) = (/ + AtS)p?{±, t) (33) 

Pp(±, t + At) = (I + A tR)p?*(±, t) (34) 
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The operator-splitting method provides the solution 
for the transport of p by making use of a Monte 
Carlo technique. In the Monte Carlo simulation, the 
density-weighted PDF at each grid cell is represented 
by an ensemble of N m stochastic elements where the 
ensemble-averaged PDF over N m delta functions re- 
places the average based on a continuous PDF. 12 


the diffusive fluxes as determined by T$/ Ax. This 
problem gets magnified if the cells also happen to be 
highly skewed. 

Such restrictions on the allowable maximum 
time step could lead to a frozen condition when the 
Monte Carlo simulation is performed with a limited 
number of stochastic particles per cell. For clarity, 
let us consider the following criterion: 


Ppm W =< PpW >= ~ W - *") ( 35 ) 

VVro n = 1 

The discrete PDF p pm (^) is defined in terms of N m 

sample values of <p n , n = 1,2, Z...N m . The statistical 

— 1/2 

error in this approximation is proportional to N m 

Using the operator-splitting method, the solu- 
tion for the PDF transport equation is obtained se- 
quentially according to the intermediate steps given 
by Eqs. (31)-(34). 

9.2 Convection /Diffusion Step 

The first step associated with convection/diffusion is 
given by: 


p£UM) = (* + A tT)pp(Tp,t) = 

<S6) 

This step is simulated by replacing a number of par- 
ticles ( = the nearest integer of c *££y m ) at <j> p {t) by 
randomly selected particles at $„(*). 

9.3 Numerical Issues Associated With Fixed 
Versus Variable Time Step 


N m > 


pAV 
c n At 


(37) 


which has to be satisfied at all grid nodes. It is es- 
timated that about 10 3 stochastic particles per cell 
are needed in order to avoid the so-called frozen con- 
dition for performing a typical 3-D gas-turbine com- 
bustor calculation. The frozen condition is referred 
to a state in which no transfer of stochastic particles 
takes place between the neighboring cells when N m 
falls below a minimum required. Scheurlen et al 27 
were the first ones to recognize the limitations asso- 
ciated with the use of a fixed time step in the Monte 
Carlo PDF computations. 

However, our experience ha s shown that this 
problem can be overcome by introducing the con- 
cept of local time-stepping which is a convergence 
improvement technique widely used in many of the 
steady-state CFD computations. In this approach, 
the solution is advanced at a variable time step for 
different grid nodes. In our present computations, it 
is determined based on 


At = min(C t fAtf , 




(38) 


It is obvious from the above equation that a nec- 
essary criterion for stability requires satisfaction of 
< 1. When the computations are performed 
with a fixed time step, this criterion tends to be too 
restrictive for most applications. Depending on the 
flow configuration, the allowable maximum time in- 
crement At is likely to be limited by a region of the 
flow field where convective fluxes dominate (such as 
close to injection holes). But in the main stream, the 
flow is usually characterized by much lower velocities. 
Resolution considerations require a higher concentra- 
tion of the grid in certain regions of the flowfield than 
in others. For example, more grid cells are clustered 
in regions where boundary layers are formed. In such 
regions the allowable maximum time increment might 
be limited in a direction dominated by the largest of 


where C%j and Ct axe calibrated constants and vrere 
assigned the values of 4 and 2.5, respectively, Atj 
is the local time step obtained from the flow (COR- 
SAIR) module, and s m \ c = The time step 

is chosen such that it permits transfer of enough par- 
ticles across the boundaries of the neighboring cells 
while ensuring that the time step used in the PDF 
computations does not deviate very much from the 
time step used in the flow solver. 

9.4 Molecular Mixing Step 

The second step associated with molecular mixing is 
given by: 

^ = -CM** ~ $*) ( 39 ) 
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The solution for this equation is updated by: 

4>t = 4>* a + (K - 4>a)e~ c * wAt (40) 

where C <j> was assigned a value of 1 . 

9.5 Spray Step 

The third step associated with the spray contribution 
is given by: 


dp<f) a 
dt 

where <f> a = h. 

The numerical solution for Eqs. (45)-(47) is in- 
tegrated by an implicit Euler scheme. 28 The result- 
ing non-linear algebraic equations are solved by the 
method of quasi-linearization. 29 


lir = ( 41 ) 

where <j> a = Va, ® — 1,2, ..., s = a — 1 

+ (42) 

where <j><j = k. The solution for the above equations 
is upgraded by a simple explicit scheme: 

.<*** _ r A At^rikmk^ 

** - e ° p AV + K (1 JKv~> 

(43) 

where a < a — 1 


9.7 Details of Combustion Chemistry 

In this section, we present an example of how com- 
bustion chemistry is handled for the case of n-heptane 
when it is modeled by a single-step global mech- 
anism of Westbrook and Dryer. 30 The correspond- 
ing rate constants in Eqs. (45)-(46) are given by 
A = 0.286 10+ 10 , a = 0.25, b = 1.25, and 

E a = 0.151 10 +O5 . This global combustion model is 
reported to provide adequate representation of tem- 
perature histories in flows not dominated by long ig- 
nition delay times. For example, the overall reaction 
representing the oxidation of the n-heptane fuel is 
given by: 

^16 + 11(02 + 3 . 76 ^ 2 )^ 


+ - 


At £ n t tn t 
pAV 


(44) 


where a = a. After a new value for enthalpy is up- 
dated, the temperature is determined iteratively from 
the solution of Eq. (10). 


7002 + 8^0 + 41 . 36^2 ( 48 ) 

Because of the constant-Schmidt-number as- 
sumption made in the PDF formulation, based on 
atomic balance of the constituent species, the mass 
fractions of N 2 , OO 2 , and H 2 0 can be shown to be 
related to the mass fractions of O 2 and CjHie by the 
following expressions: 


9.6 Reaction Step 


Vh 2 o = K 2 — K x K 2 yo 7 — K 2 yc 7 H 16 


Finally, the fourth step associated with chemical re- yco 2 = KzVh 2 o (49) 

actions is given by: 

y Ni = 1 - K 2 - K 2 K 3 - yo 2 { 1 - K X K 2 - K X K 2 K 3 )~ 


d<f*a 

dt 




where <f> a = yj. 


d<j> a 

dt 




(45) yC r H ie (1 - K 2 ~ K 2 K 3 ) 

where K x = 4.29, K 2 = 0.08943, and K 3 = 2.138. 

Using Eq. (49) results in considerable savings in 
computational time as it reduces the number of vari- 
ables in the PDF equation from five (four species and 

(46) 

one energy) to three (two species and one energy) . 


where <j> a = y a . 


9.8 Revolving Time- Weighted Averaging 
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It is noteworthy that although local time-stepping 
seems to overcome some of the problems associated 
with the PDF computations, the application of the 
Monte Carlo method requires the use of a large num- 
ber of particles because the statistical error associ- 
ated with the Monte Carlo Method is proportional 
to the inverse square root of N m , thereby making 
the use of the Monte Carlo method computationally 
very time consuming. However, a revolving averaging 
procedure used in our previous work 31 seems to alle- 
viate the need for using a large number of stochas- 
tic particles, N m) in any one given time step. Ip 
this averaging scheme, the solution provided to the 
CFD solver is based on an average of all the parti- 
cles present over the last N av time steps instead of 
an average solely based on the number of particles 
present in any one single time step. This approach 
seemed to provide smooth Monte Carlo solutions to 
the CFD solver, thereby improving the convergence 
of the coupled CFD and Monte Carlo computations. 
The reason for improvement could be attributed to 
an effective increase in the number of stochastic par- 
ticles used in the computations from N m to N av N m . 
Here, it is assumed that the solution contained within 
different iterations of the averaging procedure is sta- 
tistically independent of each other. 

10 SPRAY SOLUTION ALGORITHM 

In order to evaluate the initial conditions that are 
needed in the integration of the liquid-phase equa- 
tions, we first need to know the gas-phase properties 
at each particle location. But in order to evaluate the 
gas- phase properties, it is first necessary to identify 
the computational cell where a particle is located. It 
is a trivial task to search for the computational cell 
of the particle location in rectangular coordinates. 
However, a search for the particle location becomes 
a complicated problem when the computational cell 
is no longer rectangular in the physical domain. An 
efficient particle-search algorithm is developed and 
implemented into the Lagrangian spray solver in or- 
der to facilitate particle movement in an unstructured 
grid of mixed elements. The search is initiated in the 
form of a local search from the computational cell 
of the previous time step as the starting point. The 
location of the computational cell is determined by 
evaluating the dot product of x pc . a n = \x pc \ |a n | 
co$(<f>), where x pc is the vector defined by the dis- 
tance between the particle location and the center of 
the n-face of the computational cell, a n is the out- 
ward area normal of the n-face as shown in Fig. 2, 


and <f> is the angle between the two vectors. 



Figure 2 A vector illustration used in the particle 
search analysis. 

A simple test for the particle location requires 
that the dot product be negative over each and ev- 
ery one of the n-faces of the computational cell. If 
the test fails, the particle search is carried over to the 
adjacent cells of those faces for which the dot prod- 
uct turns out to be positive. Some of those n-faces 
might represent the boundaries of the computational 
domain while the others represent the interfaces be- 
tween two adjoining interior cells. The search is first 
carried over to the adjacent interior cells in the di- 
rection pointed out by the positive sign of the dot 
products. The boundary conditions are implemented 
only after making sure that all the possibilities lead 
to a search outside of the computational boundaries. 
This implementation ensures against any inadvertent 
application of the boundary conditions before locat- 
ing the correct interior cell. 

After the gas-phase properties at the particle lo- 
cation are known, the ordinary differential equations 
of particle position, size, and velocity are advanced by 
making use of a second-order accurate Runge-Kutta 
method. The partial differential equations govern- 
ing the droplet internal temperature distribution are 
integrated by an Euler method. Finally, after the 
liquid-phase equations are solved, the liquid-phase 
source contributions to the gas-phase equations are 
evaluated. 
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Figure 3 The flow structure of the spray code. 
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10.1 The Flow Structure of the Spray Code 
& Time- Averaging of the Interphase Source 

Terms of the Gas-Phase Equations 

The spray solver makes use of three different time 
steps: Af m / is the allowable time step, A t g i is the 
global time step, and A in is the fuel-injection time 
step. A t m i needs to be evaluated based on the small- 
est of the different time scales which are associated 
with various rate-controlling phenomena of a rapidly 
vaporizing droplet, such as those imposed by an av- 
erage droplet lifetime, the local grid spacing and a 
relaxation time scale associated with droplet veloc- 
ity among others. This restriction usually leads to 
a small time step which typically has values in the 
neighborhood of 0.01 milliseconds. However, our 
experience has shown that the convergence for the 
steady-state computations could be improved greatly 
by supplying the flow and PDF solvers with the in- 
terphase terms obtained from a time-averaging proce- 
dure, where the averaging is performed over an aver- 
age lifetime of the droplets, A t g i. The variable, A t g ^ 
has values in the neighborhood of 1 ms. 

The averaging scheme could be explained better 
through the use of a flow chart shown in Fig. 3. The 
main spray solver is invoked by a call to the control- 
ling routine which executes the following steps: 

1. It first initializes the source terms to zero. 

2. Checks to see if new particles need to be 
introduced. 

3. Advances liquid-phase equations over a pre- 
specified time step, A t m j, with calls to the 
following routines: 

- Does a particle search and assigns par- 
ticles based on the parallel strategy im- 
plemented. 

- Interpolates gas-phase properties at 
the particle location. 

- Advances liquid-phase equations and, 
also, deletes any particles that are no 
longer needed in the computations. 

4. Evaluates the liquid-phase source-term con- 
tributions, Smh f° r use i n the gas-phase 
equations. 

5. Continues with steps (2) and (3) until the 
computations are completed over a global 
time step of A t g j. 


6. Returns control to other solvers, e.g., flow 
or PDF, and supplies them with source 
terms, S g i, averaged over At g i. 

The time-averaged contribution of these source 
terms, S g i> is given by: 

M A 

s s > = £ tfr s- ( 5 °) 

m=l 

where 

u 

A< mi = A tgi (51) 

m— 1 


11 COUPLING BETWEEN THE THREE 
SOLVERS 

For the PDF solver, the mean gas-phase velocity, tur- 
bulence diffusivity and frequency are provided as in- 
puts from the CFD solver and the modeled spray 
source terms from the liquid-phase solver. And, in 
turn, the Monte-Carlo solver supplies the tempera- 
ture and species fields to the other two solvers. The 
CFD code also receives the liquid-phase source terms 
as inputs from the spray solver. For the spray solver, 
the needed gas-phase velocity and scalar fields are 
supplied by the other two solvers. The liquid-phase, 
PDF, and CFD solvers are advanced sequentially in 
an iterative manner until a converged solution is ob- 
tained. It should also be noted that both the PDF 
and spray solvers are called once at every specified 
number of CFD iterations. All three PDF, CFD, and 
spray codes were coupled and parallelized in such a 
way in order to achieve maximum efficiency. 

12 RESULTS AND DISCUSSION 

The calculation procedure was applied to pre- 
dict the flow properties of an unconfined spray flame 
for which some experimental data was reported by 
McDonnel and Samuelson. 20 We also discuss the re- 
sults from a 3D test case which was designed to 
demonstrate the viability of the present solution pro- 
cedure for use in 3D computations. 

12.1 Unconfined Spray Flame 
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Fig. 7 Global structure of a spray flame. 
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Fig. 9 Oxygen mass fraction contours. 
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A detailed experimental data set for the case of 
an unconfined methanol spray was reported by Mc- 
Donnel and Samuelson. 20 The schematic of the ex- 
perimental facility used at UCI (University of Cali- 
fornia, Irvine) is shown in Fig. 4. It made use of a 
Parker and Hannifin RSA (Research Simplex Atom- 
izer) whose simplex tip has a flow number of 1.36 and 
a nominal spray angle of 35 degrees. The reported 
methanol and air mass flow rates were 1.26 and 1.32 
g/s, respectively. The spray was injected downwards 
from the center of a 495 x 495 mm square duct and 
air was pulled through the top of the duct by a blower 
at a bulk velocity of 0.8 m/s in order to provide ade- 
quate entrainment needs. Both the droplet and gas- 
phase velocities as well as the droplet sizes were mea- 
sured by making use of a two-component PDI (Phase 
Doppler Interferometry), and the gas-phase temper- 
atures were measured by using a traversing hot-wire 
thermocouple. Using the setup shown in Fig. 4, sev- 
eral measurements involving the gas-phase velocity, 
droplet size and velocity, droplet number flux, and 
mean gas-phase temperatures were reported at dif- 
ferent axial locations starting from 2.5 cm. 

A photograph of the burning spray is shown in 
Fig. 5. For the case considered, it produces a narrow 
illuminating flame at the center. 

For the axisymmetric configuration, both the 
PDF and non-PDF computations were performed on 
a 2D grid of 1850 triangular elements as shown in 
Fig. 6. Since in complex 3D geometries, it is not al- 
ways possible to identify and delineate critical regions 
of a flow-field, a relatively coarse mesh was chosen to 
see how well the flow-field could be computed without 
resorting to a fine grid. The turbulent Schmidt and 
Prandtl numbers were taken to have a value of 0.70. 
The PDF solution was obtained by making use of 100 
particles per cell. The temperature and species fields 
supplied to the CFD and liquid-phase solvers were 
obtained from averaging the PDF solutions over a pe- 
riod of the last 100 time-steps. The calculations were 
advanced until a steady state solution was reached 
by making use of the following time steps: A t g was 
determined based on a local time-step with a CFL 
number of 1, At injection = 1-0 ms, and A tk = 0.01 
ms. At the end of every liquid-phase injection time 
step, a new droplets group was introduced. 

12.1.1 Global Features of the Spray Flame 

The global features of the spray flame are shown 
in Figs. 7a and 7b by presenting a composite view of 


the droplet locations and the mean gas-phase temper- 
atures and velocity vectors. The filled white circles 
show the location of the droplets. The droplet sizes 
range from few microns to 140 microns. The shaded 
contour lines show the temperature distribution, and 
the arrows denote the velocity vectors. Fig. 7a shows 
the results from the PDF method and Fig. 7b from 
the non-PDF (conventional CFD) method. 

First, let us first look at the droplet distribution. 
As expected, because of the prevailing high temper- 
atures, the droplets in the central region of the spray 
tend to vaporize faster than those present in the outer 
regions of the spray. For that reason, the average size 
of the droplets in the central region is much lower 
than those present elsewhere. However, most of the 
droplet mass is contained within the droplets of the 
high temperature region. The largest droplets found 
outside of the high-temperation region are sometimes 
called in the literature as rogue droplets. 20 

Next, let us look at the results from the PDF 
computations. Combustion seems to be initiated by 
a flame front stabilized in the lower velocity region of 
the outer shear layer. Starting from there, the high 
temperature region spreads in a long v-neck shape 
as a result of two distinct flames being formed. The 
inner flame has the characteristics of a turbulent pre- 
mixed flame and the outer flame shows the character- 
istics of a turbulent diffusion flame. The distinction 
between the two flames becomes more obvious when 
we will look at the fuel and oxygen mass-fraction con- 
tour plots shown in Figs. 8 and 9. Methanol is known 
to vaporize rapidly, its vapor has the same density as 
air, and its liquid saturation temperature is approx- 
imately around 263 K which can easily be achieved 
prior to atomization. 20 Early vaporization and mix- 
ing of fuel vapor with air leads to the formation of an 
inner premixed flame. Further vaporization leads to 
a large accumulation of fuel vapor inside of the cen- 
tral high-temperature region (Fig. 8a). This region 
is also devoid of oxygen (Fig. 9a). Here, combustion 
takes place with the formation of a diffusion flame 
where the fuel from the inner central region mixes 
and burns with the surrounding air from the outer 
region. However, it is noteworthy that in sprays this 
kind of characterization of attributing the flames to 
be either premixed or diffusion could only be made 
in a general sense. Because of vaporizing droplets 
present within an active flame zone, it makes any 
kind of strict characterization meaningless. 

On the other hand, the non-PDF computations 
in Fig. 7b show for the combustion to occur in a pre- 
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Monte Carlo particle temperatue, K Monte Carlo particle temperatue, 





dominantly vaporization-controlled reaction regime. 
As a result, the high temperature region is spread 
over a wider region. It lacks a well-defined flame 
structure that was observed with the PDF compu- 
tations. 

Near the centerline, the PDF results show for 
the jet to retain its axial momentum better than the 
non-PDF solution. The PDF results also show less 
radial spreading of the jet further downstream. As we 
sill see later, the PDF results are more in compliance 
with the experimental data than those predicted by 
the non-PDF calculations. 

12.1.2 Gas-Phase Mass Fraction Contours 

Figs. 8 an 9 show the mass fraction contours of 
methanol and oxygen, respectively. First, looking at 
the methanol mass fractions reinforces the remarks 
that were made earlier on the PDF results. Early 
vaporization leads to an accumulation of fuel vapor 
in the inner core region of the jet which supports a 
premixed inner flame (Fig. 7a). But, further down- 
stream, there is a large accumulation of fuel vapor 
in the the region where high-temperature products 
are present. This fuel-rich region is also devoid of 
oxygen (Fig. 9a). Mixing of this fuel vapor with 
the surrounding outer air supports a diffusion flame. 
However, this outer flame also feeds on the fuel va- 
por from the vaporizing droplets present in its active 
combustion region. The oxygen mass fractions (Fig. 
9a) also show clearly a well-defined flame structure. 

However, the non-PDF computations show a 
slight accumulation of fuel vapor in only a small re- 
gion of the entine domain which implies instanta- 
neous burning of the vaporized fuel. Unlike the PDF 
predictions, the oxygen mass fractions lack a well- 
defined flame structure as combustion seems to occur 
over a much wider region in a vaporization-controlled 
regime. 

12.1.3 Gas-Phase Monte Carlo Particle Temperature 
Distribution 

In order to examine further what is causing the 
major differences found in the PDF and non-PDF 
results, let us take a look at the scatter plots of the 
Monte Carlo particle temperature distribution shown 
in Figs. lOa-b. Fig. 10a contains the particle dis- 
tribution at four different radial nodes of an axial 
location at 0.06 m, and Fig, 10b provides similar in- 
formation at an axial location of 0.12 m. The particle 
temperatures shown here have been pre-sorted to fall 
in a descending order. 


Referring to Fig. 7a, it is obvious that the four 
nodes at the first axial location are chosen to examine 
the differences in the Monte Carlo particle distribu- 
tion associated with the inner and outer flames. The 
second axial location is chosen to examine the differ- 
ences further downstream. 

The first node is from the central region region 
of the jet, and it falls in the the outer (reactants) 
region of the inner (premixed) flame. The second and 
third nodes are located close to the active flame zones 
of the inner and outer flames, respectively. And the 
fourth one falls in the outer (surrounding air) region 
of the outer (diffusion) flame. 

In the first node, few of the particles are at flame 
temperature and most of the rest are at a lower tem- 
perature of the surrounding reactants. In the next 
cell, most of the particles are at a flame temperature 
and it is followed by an abrupt drop near the end. 
The abrupt changes in the temperature distribution 
at the first two locations are indicative of the inter- 
mittency effects associated with a turbulent premixed 
flame. It is obvious that the inner flame is character- 
ized by large fluctuations in temperature. If the cor- 
responding PDFs were constructed, the PDFs would 
show two distinct peaks of varying amplitudes de- 
pending upon the proximity to the flame front. One 
of the peaks is centered around the flame temperature 
and the other at the temperature of the surrounding 
reactants. 

The particles distribution in the next two cells 
are associated with the outer diffusion flame. Unlike 
the inner flame, the particle distribution for the outer 
flame is mainly characterized by a single peak with 
the temperatures ranging between 2100 to 1400 K at 
the third cell and 1000 to 900 K at the fourth cell. 
However, the PDF at the third location (near the 
flame) would show a distribution over a wider tem- 
perature range with a relatively large variance com- 
pared to one associated with the fourth location. 

Figure 10b shows a small variance at the first 
two cells with a mean closer to the flame tempera- 
ture. The small variance in the temperature is due to 
the changes in the composition of the product gases 
which are caused by the addition of fuel vapor from 
the vaporizing droplets. The particle distribution in 
the next two nodes is essentially what one expects 
to see in the outer surrounding regions of a typical 
diffusion flame. As expected, the mean temperature 
shows the inverse functional relationship with the in- 
creasing normal distance away from the flame front. 
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Fig. 1 1 Gas-phase mean temperature comparisons. 
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Fig. 13 Gas-phase mean radial velocity comparisons (experimental rms is shown as an error bar). 
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12.1.4 Mean Gas-Phase Temperature Comparisons 

Figs, lla-d show the comparisons for the radial 
profiles of mean gas temperature at four different ax- 
ial locations. Both the PDF and non-PDF computa- 
tions seem to predict higher peak temperatures than 
those measured. Some of the differences found could 
be attributed for the following reasons: 

(1) The wetting of the thermo-couples might 
have contributed to some of the uncertaininty ob- 
served in the the reported temperatures. For this 
reason, the comparisons in this section are mainly 
meant to provide a qualitative description. 

(2) The uncertainty contained in the experimen- 
tal data is not clear because no error bands were pro- 
vided for the data reported. 

(3) Even though ethanol flame temperature is 
about 2100 K under normal conditions, the maxi- 
mum attainable temperature is also a function of the 
initial temperature of the gas mixture. The bulk in- 
flow temperature of about 500 K is higher than the 
normal ambient temperature which leads to a corre- 
spondingly higher flame temperature. 

(4) The use of a single-step global mechanism 
is known to overpredict the flame temperatures by 
about 100 to 200 K. 

(5) Some of the temperature drop could be at- 
tributed to the higher emissivity and radiation cool- 
ing rate of the spray combustion flame. 

At the inflow in Fig. 11a, if one is wonder- 
ing why should there exist any differences at all be- 
tween the measured and computed (both PDF and 
non-PDF) temperatures, it is because the predictions 
represent the temperatures extrapolated from those 
known at the cell centers instead of those specified 
on the boundary faces of the cells where the inflow 
conditions were applied. In Fig. 11a, the PDF re- 
sults show very little change from the specified inflow 
temperatures, and the results are very close to the 
specified the experimental data. However, the non- 
PDF computations show an abrupt jump from the 
specified inflow temperatures. 

At the next location in Fig. lib, the non- 
PDF shows the formation of a hot central region with 
the centerline peak reaching a flame temperature of 
about 2200 K. Unlike the non-PDF, both the mea- 
sured and the PDF show for the temperature to peak 
in the outer region of the spray, and the peak tem- 
perature location is correctly predicted by the PDF 
computations. However, the PDF underpredicts the 
centerline temperature while overpredicting the peak 


temperature. If you recall from our earlier discus- 
sion on the temperature fluctuations in this region, 
it is quite likely that the measuring devices (thermo- 
couples) may have difficulty in capturing the correct 
mean temperatures in this region. Farther away from 
the outer region of the spray, there is no physical ex- 
planation as to why the measured temperatures are 
higher than those predicted by the PDF other than 
any differences stemming from the neglect of radia- 
tive heat transfer. 

At the third location, near the centerline, the 
PDF results are in compliance with the measured 
data. At the last two locations, the PDF again 
correctly predicts the location of the peak temper- 
atures. But, in both locations, the PDF clearly over- 
predicts the peak temperatures while underpredicting 
the temperatures in the outer regions of the spray. 
On the other hand, the non-PDF results show a fur- 
ther broadening of the central high-temperature re- 
gion and a greater radial temperature spreading into 
the outer region. 

12.1.5 Gas- Phase Velocity Comparisons 

Figs. 12a-d show the comparisons for the mean 
gas-phase axial velocities at four different axial lo- 
cations. As expected near the inflow, both the pre- 
dictions and measurements show similar behavior as 
shown in Fig. 12a. But at the next three downstream 
locations of Figs. 12b-d, both the PDF and non- 
PDF predictions underpredict the velocities near the 
centerline. The PDF predictions are more closer to 
what was observed experimentally, and the compar- 
isons become progressively worse further downstream 
with the non-PDF computations. 

The radial velocity comparisons are shown Figs. 
13a-d. It is noteworthy that the experimental data 
seem to show a great deal of noise when it comes to 
measuring the radial velocities. For that reason, the 
rms component of the radial velocity is superimposed 
as a vertical error bar. The error bars clearly show 
that the fluctuations are indeed very large and in 
some instances even exceed the corresponding mean. 
Similar noise was reported in the data used in our pre- 
vious comparisons of a swirl-stabilized spray case. 9,32 
The PDF correctly predicts the locations of the peaks 
at all four locations but seem to underpredict the 
magnitudes by a considerable measure. However, the 
non-PDF computations tend to overpredict the radial 
distances of the peak, and the comparisons become 
progressively worse further downstream. 
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Fig. 15 Drop mean radial velocity comparions (experimental rms is shown as an error bar). 
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Fig. 16 Scatter plot of drop size comparisons. 





Fig. 17 Schematic of a 3D test case. 
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Fig. 19 Temperature distribution of a 3D test case. 


29 















12.1.6 Droplet Velocity fc Size Comparisons 

The scatter plots of the mean droplet axial ve- 
locity are presented in Figs. 14a-d at four different 
axial locations. These plots include all of the reported 
measured velocities as well as the Lagrangian par- 
ticle velocities. The experimental reported data is 
only limited to measurements in the radial direction 
of 0.02 m from the centerline distance. 

The PDF seems to provide a better agree- 
ment than the non-PDF, and the non-PDF com- 
parisons become progressively worse further down- 
stream. These comparisons are consistent with the 
mean gas-phase velocity comparisons discussed ear- 
lier. The differences observed near the centerline with 
the PDF predictions could partially be attributed to 
the corresponding behavior observed earlier in the 
gas-phase velocities. As we have discussed earlier, 
because of the high temperatures that existed in the 
central region of the jet, the droplet sizes on the av- 
erage in that region tend to be smaller because they 
vaporize faster. This would in turn cause the droplets 
to relax toward the surrounding gas-phase velocities 
much faster because of the drag forces acting on the 
smaller droplets. 

Figs. 15a-d show the scatter plots of mean 
droplet radial velocity. As in the gas-phase, the ex- 
perimental data show a great deal of noise when it 
comes to measuring the radial velocities. The error 
bars based on the rms component clearly show that 
the fluctuations are very large and in some instances 
even exceed their corresponding mean. Both the PDF 
and non-PDF comparisons underpredict the exper- 
imental data. The PDF computations are able to 
capture the qualitative trends correctly but the non- 
PDF comparisons become progressively worse further 
downstream. 

Fig. 16a-d show the scatter plots of droplet 
sizes. The experimental data is reported in the form 
of several radial measurements for a given droplet- 
size range. Unfortunately, the plots do not contain 
any information on the number density of a given par- 
ticle group. Although the experimental data shows 
a wider presence radially, only some of the droplet 
groups do contain most of of the mass. The results 
show that the droplet sizes are well represented by 
both the PDF and non-PDF computations but the 
non-PDF computations show a far less number of 
droplets within the reported experimental range of 
the last location. In the non-PDF computations, the 
droplets traversing the central region of the jet vapor- 


ize more rapidly as they are exposed to temperatures 
of about 2100 K for a longer period of time. This 
would explain the reason for the lack of enough num- 
ber droplets in this range. 

12,2 A 3D Test Case 

Here, the results of a 3D test case are summa- 
rized. Fig. 17 shows the geometry of the combustor 
configuration used in the computations. The inflow 
air flows through an annular opening of 4 5- degree 
swirl vanes with a flow speed of about 28.25 m/s and 
a mass flow rate of about 41.1 g/s. For the annular 
opening, the dimensions of the outer and inner radii 
are 3.2 and 1.5 cm, respectively. N-heptane fuel was 
injected in the form of single hollow cone stream with 
a mass flow rate of about 2.7 g/s. Solid wall bound- 
ary conditions were applied on the boundaries shown 
with the surface grids; the top and bottom surfaces 
as well as those on the frontal surface. The flow exits 
through the backward surface and periodic conditions 
were applied on the remaining left and right surfaces. 
Both the PDF and the non- PDF computations were 
performed for this case on a computational grid com- 
prising of 8430 tetrahedral elements and 100 Monte 
Carlo PDF particles per cell. 

The velocity vector plot of Fig. 18 shows the 
formation of a vortex swirling in the clock-wise di- 
rection near the inflow boundary. The flow becomes 
more unidirectional further downstream as it flows to- 
wards the exit boundary. Initially, most of the spray 
particles move towards the upper left corner near the 
top wall, and from there onwards, the particles seem 
to drift downstream towards the exit. This leads to a 
concentration of fuel vapor near that left location of 
the upper wall from the vaporizing droplets as shown 
in Fig. 20. The temperature contours of Fig. 19 
show that a flame originating near the left corner of 
the top wall surface, propagates radially outwards as 
it moves downstream. The fuel and oxygen mass frac- 
tions of Figs. 20 and 21 confirm this observation as 
the fuel vapor from the vaporizing droplets near the 
upper wall mixes and burns with surrounding air to 
support a flame originating from that location. 

13 CONCLUDING REMARKS 

A solution procedure was developed for spray 
computations based on the scalar Monte Carlo PDF 
method for application with unstructured grids and 
parallel computing. The solution procedure seems 
to be able to capture the overall structure of a 
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Table 1 

. Cpu time (sec) per cycle versus number of PEs. 



Number of processors 

Solver 

Characteristic 

2 

5 

10 

CFD 

5 steps/cycle 

2.50 

1.25 

0.75 

PDF 

1 step/cycle 

6.5 

2.9 

1.9 

Spray 

100 steps/cycle 

1.70 

0.64 

0.53 


spray flame well, and its application to several 
spray flames (both confined as well as unconfined) 
showed reasonable agreement with the available spray 
measurements . 7> 9 

The detailed comparisons made for the case of 
a reacting spray illustrated the importance of chem- 
istry/turbulence interactions in the modeling of a re- 
acting spray. The PDF results were found to be closer 
to the reported experimental data when compared 
with the non-PDF solution. The PDF computa- 
tions predict that some of the combustion occurs in a 
predominantly premixed-flame environment and the 
rest in a predominantly diffusion-flame environment. 
However, the non-PDF solution predicts wrongly for 
the combustion to occur in a vaporization-controlled 
regime. Near the premixed flame, the Monte Carlo 
particle temperature distribution shows two distinct 
peaks: one centered around the flame temperature 
and the other around the surrounding-gas tempera- 
ture. Near the diffusion flame, the Monte Carlo par- 
ticle temperature distribution shows a single peak. In 
both cases, the PDF’s shape and strength are found 
to vary substantially depending upon the proximity 
to the flame surface. The results cast some ambiq- 
uity regarding the applicability of the widely used 
assumed-shape PDF methods in spray computations. 

The paper also demonstrated the results from a 
3D test case which wets designed primarily to demon- 
strate the computational viability of the Monte Carlo 
PDF method for 3D computations and also to ex- 
amine the applicability of the spray particle search 
algorithm in 3D computations and of the newly- 
implemented periodic boundary conditions. 

As mentioned in the previous section, the com- 
putations for the 3D test case were performed on a 
grid of 8,450 elements with a total of 0.845 million 
particles (=100 particles/cell). The computations 
were performed on one of the NASA Ames Research 
Center’s parallel computer platforms called Turing 


which is a SGI Origin work-station with 24 PEs (Pro- 
cessor Elements). Table 1 summarizes the cpu times 
per cycle taken by the PDF, spray, and CFD solvers 
vs the number of PEs. All of the CFD, PDF and 
spray solvers show good parallel performance with 
an increase in the number of processors. 

The additional computational burden associ- 
ated with use of a Monte Carlo PDF method should 
be considered to be moderate in view of the the 
normal cpu times taken by the conventional single- 
processor CFD/spray solvers. The viability of the 
present method for its application to the model- 
ing of practical combustion devices is demonstrated 
through the use of unstructured grids and the ability 
to run the computations on massively parallel com- 
puters. 
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